#include <cmath>
#include <iostream>
#include "StencilMath.h"
#include "swt_env_win_linux.h"
#include "lct_env_macros.h"

#ifdef _SWT_WIN32_
# define M_PI 3.14159265358979323846
#endif


#define INDEX(x, y, z)  ((x) + (__int64) ldx * ((y) + (__int64) ldy * (z)))


void StencilMath::stencil_order_4_kern(double* dst, const double* src, const double* coeffc2_dist,
                                       int start_x, int start_y, int start_z,
                                       int end_x, int end_y, int end_z, int ldx, int ldy) {
    int i, j, k;

#ifdef _OMP_H_
    #pragma omp parallel for schedule(dynamic)
#endif
    for(k = start_z; k < end_z; k++) {
        for(j = start_y; j < end_y; j++) {

            for(i = start_x; i < end_x; i++) {
                double coeffc2 = coeffc2_dist[INDEX(i, j, k)];

                dst[INDEX(i, j, k)] = -dst[INDEX(i, j, k)]
                                      + (2 - 7.5 * coeffc2) * src[INDEX(i, j, k)]
                                      + coeffc2 * (src[INDEX(i - 1, j, k)] * 4 / 3
                                                   + src[INDEX(i + 1, j, k)] * 4 / 3
                                                   - src[INDEX(i - 2, j, k)] / 12
                                                   - src[INDEX(i + 2, j, k)] / 12

                                                   + src[INDEX(i, j - 1, k)] * 4 / 3
                                                   + src[INDEX(i, j + 1, k)] * 4 / 3
                                                   - src[INDEX(i, j - 2, k)] / 12
                                                   - src[INDEX(i, j + 2, k)] / 12

                                                   + src[INDEX(i, j, k - 1)] * 4 / 3
                                                   + src[INDEX(i, j, k + 1)] * 4 / 3
                                                   - src[INDEX(i, j, k - 2)] / 12
                                                   - src[INDEX(i, j, k + 2)] / 12
                );
            }
        }
    }
}


void StencilMath::stencil_order_4_pml_kern(double* dst, const double* src, const double* coeffc2_dist, const double* c2dtQ_dist,
                                           int start_x, int start_y, int start_z,
                                           int end_x, int end_y, int end_z, int ldx, int ldy) {
    int i, j, k;

#ifdef _OMP_H_
    #pragma omp parallel for schedule(dynamic)
#endif
    for(k = start_z; k < end_z; k++) {
        for(j = start_y; j < end_y; j++) {

            for(i = start_x; i < end_x; i++) {
                double Qc2dt = c2dtQ_dist[INDEX(i, j, k)];
                double coeffc2 = coeffc2_dist[INDEX(i, j, k)];

                dst[INDEX(i, j, k)] = ((Qc2dt - 1) * dst[INDEX(i, j, k)]
                                       + ((2 - Qc2dt * Qc2dt) - 7.5 * coeffc2) * src[INDEX(i, j, k)]
                                       + coeffc2 * (src[INDEX(i - 1, j, k)] * 4 / 3         // line 1
                                                    + src[INDEX(i + 1, j, k)] * 4 / 3       // line 2
                                                    - src[INDEX(i - 2, j, k)] / 12          // line 3
                                                    - src[INDEX(i + 2, j, k)] / 12          // line 4

                                                    + src[INDEX(i, j - 1, k)] * 4 / 3       // line 5
                                                    + src[INDEX(i, j + 1, k)] * 4 / 3       // line 6
                                                    - src[INDEX(i, j - 2, k)] / 12          // line 7
                                                    - src[INDEX(i, j + 2, k)] / 12          // line 8

                                                    + src[INDEX(i, j, k - 1)] * 4 / 3       // line 9
                                                    + src[INDEX(i, j, k + 1)] * 4 / 3       // line 10
                                                    - src[INDEX(i, j, k - 2)] / 12          // line 11
                                                    - src[INDEX(i, j, k + 2)] / 12          // line 12
                                                   )
                                      ) / (1 + Qc2dt);
            }


        }
    }
}


double StencilMath::Ricker(double f, double t, double t0) {
    double phi = M_PI * f * t - 4. - t0;     // -4 works as t0 to smoothen the amplitude curve
    double result = (1 - 2 * phi * phi) * exp(-phi * phi);

    return result;
}
